home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / fits.i < prev    next >
Text File  |  1996-02-21  |  26KB  |  878 lines
  1. /* 
  2.  * fits.i --
  3.  *
  4.  *    Routines to manipulate FITS files (IAU astronomical data format).
  5.  *    Provides functions:
  6.  *      - fitsAddComment: add comments in a FITS header;
  7.  *      - fitsAddHistory: add history records in a FITS header;
  8.  *      - fitsBitpix:     compute bit per pixel;
  9.  *      - fitsFixHeader:  check/fix FITS header;
  10.  *      - fitsHeader:     provides default FITS header,
  11.  *      - fitsRead:       read FITS file;
  12.  *      - fitsRescale:    rescale data into integer array;
  13.  *      - fitsWrite:      write FITS file.
  14.  *
  15.  * $Id: fits.i,v 1.5 1996/02/21 10:55:30 eric Exp $
  16.  *
  17.  * Copyright (c) 1996, Eric THIEBAUT (thiebaut@obs.univ-lyon1.fr, Centre de
  18.  * Recherche Astrophysique de Lyon, 9 avenue Charles  Andre,  F-69561 Saint
  19.  * Genis Laval Cedex).
  20.  *
  21.  * This program is free software; you can redistribute it and/or  modify it
  22.  * under the terms of the GNU General Public License  as  published  by the
  23.  * Free Software Foundation; either version 2 of the License,  or  (at your
  24.  * option) any later version.
  25.  *
  26.  * This program is distributed in the hope  that  it  will  be  useful, but
  27.  * WITHOUT  ANY   WARRANTY;   without   even   the   implied   warranty  of
  28.  * MERCHANTABILITY or  FITNESS  FOR  A  PARTICULAR  PURPOSE.   See  the GNU
  29.  * General Public License for more details (to receive a  copy  of  the GNU
  30.  * General Public License, write to the Free Software Foundation, Inc., 675
  31.  * Mass Ave, Cambridge, MA 02139, USA).
  32.  */
  33.  
  34. /*
  35.  * History:
  36.  *    Known bugs in release 1.1:
  37.  *      - sometimes fails to compile with Yorick-1.1a on a HP workstation;
  38.  *      - bugs in `fitsWrite'.
  39.  *
  40.  *    02/16/96 Eric THIEBAUT:
  41.  *      - fixed bugs in `fitsWrite' when array is of integer type;
  42.  *      - added some functions and cleanup the code;
  43.  *      - removed spurious debug messages;
  44.  *      - seems to properly compile with Yorick-1.2.
  45.  *
  46.  *    02/16/96 release 1.3.
  47.  *
  48.  *    02/17/96 Eric THIEBAUT:
  49.  *      - when rescaling in fitsRead, bitpix=32 -> double, otherwise float;
  50.  *      - fitsRead: set file primitive to be that of a SUN (most significant
  51.  *        byte first) to follow FITS standard and make data portable (that
  52.  *        was already done in fitsWrite and was fixed by Dave in its 1.2
  53.  *        release...);
  54.  *      - fitsRead: by default ignores zero *AND* unit dimensions, this
  55.  *        was not what the documentation pretends; fix it to follow
  56.  *        documentation;
  57.  *      - incorporates most of the changes made by David Munro.
  58.  *
  59.  *    02/17/96 release 1.4.
  60.  *
  61.  *    02/21/96 Eric THIEBAUT:
  62.  *      - in fitsRead: speed-up parsing of FITS header (a side effect is
  63.  *        that the maximum number of dimensions is 33);
  64.  *      - in fitsWrite and _fitsPut...: use extern variables to simplify
  65.  *        (and speed-up?) the writting of the header;
  66.  *
  67.  *    02/21/96 release 1.5.
  68.  *
  69.  *    Known bugs in release 1.5:
  70.  *      - Even if fits_max_naxis is changed, the maximum number of
  71.  *        dimensions is 33.  I do not think that this is a limitation
  72.  *        thus I won't fix this ;-).
  73.  */
  74.  
  75. require, "string.i";
  76.  
  77. /*----------------------------------------------------------------------*/
  78.  
  79. local fits_max_naxis;
  80. /* DOCUMENT fits_max_naxis -- maximum number of allowed axis in FITS
  81.     Its value should *NOT* be modified from a Yorick session.
  82. */
  83. fits_max_naxis= 9;    /* if you need more axis, change this value in this
  84.              * file and source it again. */
  85.  
  86. struct FitsHeader {
  87.   int        bitpix;        //    8: pixels are unsigned bytes
  88.                 //   16: pixels are signed 2-byte integers
  89.                 //   32: pixels are signed 4-byte integers
  90.                 //  -32: pixels are 4-byte floating points
  91.                 //  -64: pixels are 8-byte floating points
  92.                 // -128: pixels are 2x8-byte complex
  93.   int        naxis;        // number of axis
  94.   int        axis(fits_max_naxis);    // number of pixel along axis No. n
  95.   // Pixel intensity / brighness:
  96.   double    bscale;        // pixelValue = BZERO+BSCALE*fileValue
  97.   double    bzero;        // pixelValue = BZERO+BSCALE*fileValue
  98.   string    bunit;        // brightness unit
  99.   double    datamax;    // maximum data value in the file
  100.   double    datamin;    // minimum data value in the file
  101.   // Miscellaneous information:
  102.   string    object;        // image name
  103.   string    date;        // date of file creation (dd/mm/yy)
  104.   string    date_obs;    // date of data acquisition (dd/mm/yy)
  105.   string    origin;        // institution
  106.   string    instrume;    // data acquisition instrument
  107.   string    telescop;    // data acquisition telescope
  108.   string    observer;    // observer name/identification
  109.   string    history;    // newline separated history lines
  110.   string    comment;    // newline separated comment lines
  111.   // coordinate system for each axis:
  112.   double    epoch;        // epoch of coordinate system (year)
  113.   double    crval(fits_max_naxis);    // coord = CRVAL+(pixel-CRPIX)*CDELT
  114.   double    crpix(fits_max_naxis);    // coord = CRVAL+(pixel-CRPIX)*CDELT
  115.   double    cdelt(fits_max_naxis);    // coord = CRVAL+(pixel-CRPIX)*CDELT
  116.   string    ctype(fits_max_naxis);    // type of physical coordinate
  117.   double    crota(fits_max_naxis);    // rotation angle of axis No. #
  118. }
  119.  
  120. /*----------------------------------------------------------------------*/
  121.  
  122. func fitsHeader(&header)
  123. /* DOCUMENT header= fitsHeader();
  124.             fitsHeader, header;
  125.     Creates  a  structure  FitsHeader  with  defaults.   Equivalent to:
  126.         header= FitsHeader(bscale= 1., bzero= 0.);
  127.  
  128.    SEE ALSO: fitsRead, fitsWrite.
  129. */
  130. {
  131.   return (header= FitsHeader(bscale= 1., bzero= 0.));
  132. }
  133.  
  134. func fitsFixHeader(&header)
  135. /* DOCUMENT fitsFixHeader, header
  136.     Check consistency of FITS header.
  137.  
  138.    SEE ALSO: fitsRead, fitsWrite.
  139. */
  140. {
  141.   if (structof(header) != FitsHeader)
  142.     error, "bad type for HEADER (not FitsHeader)";
  143.   if (header.datamin > header.datamax)
  144.     error, "DATAMIN greater than DATAMAX in FITS header!";
  145.   if (header.bscale == 0. && header.bzero != 0.) {
  146.     write, "WARNING: BZERO ignored because BSCALE is zero";
  147.     header.bzero= 0.;
  148.   }
  149. }
  150.  
  151. func fitsAddComment(&header, str)
  152. /* DOCUMENT fitsAddComment, header, string
  153.     Add STRING as a comment in HEADER  (a  FitsHeader  structure).   If
  154.     STRING is too long, it is broken in pieces.
  155.  
  156.    SEE ALSO: fitsRead, fitsWrite.
  157. */
  158. {
  159.   str= strtrim(str);    // remove trailing/leading blanks
  160.   while (strlen(str)) {
  161.     header.comment += strpart(str, :69)+"\n";
  162.     str= strtrim(strpart(str, 70:), 1);
  163.   }
  164. }
  165.  
  166. func fitsAddHistory(&header, str, stamp=)
  167. /* DOCUMENT fitsAddHistory, header, string
  168.     Add STRING as an history card in HEADER  (a  FitsHeader structure).
  169.     If STRING is too long, it is broken in pieces.   Keyword  STAMP may
  170.     be used with a  zero  value  to  avoid  the  additional  time stamp
  171.     (default is to add a time stamp).
  172.  
  173.    SEE ALSO: fitsRead, fitsWrite.
  174. */
  175. {
  176.   str= strtrim(str);    // remove trailing/leading blanks
  177.   if (is_void(stamp))
  178.     stamp= 1N;    // default is to add a time stamp
  179.   if (stamp) {
  180.     str= getdate()+" "+gettime()+": "+str;
  181.     indent= "                   ";
  182.   } else {
  183.     indent= "";
  184.   }
  185.   len= 69 - strlen(indent);
  186.   header.history += strpart(str, :69)+"\n";
  187.   str= strtrim(strpart(str, 70:), 1);
  188.   while (strlen(str)) {
  189.     header.history += indent + strpart(str, :len) + "\n";
  190.     str= strtrim(strpart(str, len+1:), 1);
  191.   }
  192. }
  193.  
  194. func fitsBitpix(data)
  195. /* DOCUMENT fitsBitpix(data)
  196.     Return bits per pixel for array data.
  197.    
  198.    SEE ALSO: fitsWrite.
  199. */
  200. {
  201.   s= structof(data);
  202.   if (s == char)
  203.     return 8 * sizeof(char);
  204.   if (s == short)
  205.     return 8 * sizeof(short);
  206.   if (s == int)
  207.     return 8 * sizeof(int);
  208.   if (s == long)
  209.     return 8 * sizeof(long);
  210.   if (s == float)
  211.     return -8 * sizeof(float);
  212.   if (s == double)
  213.     return -8 * sizeof(double);
  214.   if (s == complex) {
  215.     write, "WARNING: using BITPIX=-128 is not standard FITS";
  216.     return -8 * sizeof(complex);
  217.   }
  218.   error, "bad data type for DATA";
  219. }
  220.  
  221. /*----------------------------------------------------------------------*/
  222.  
  223. func fitsRescale(data, bitpix, &bscale, &bzero, data_min=, data_max=)
  224. /* DOCUMENT rescaled_data= fitsRescale(data, bitpix, bscale, bzero)
  225.  
  226.     Linearly rescale  the  values  of  input  array  DATA  to  fit into
  227.     integers with BITPIX  bits  per  value  (i.e.,  `char',  `short' or
  228.     `long' for BITPIX being 8, 16 and 32 respectively).
  229.  
  230.     Arguments BSCALE and BZERO are optional and  purely  outputs passed
  231.     by address.  Their value will be set so  that:
  232.         DATA(i) = BZERO + BSCALE * RESCALED_DATA(i)
  233.     where the equality may not be exact due to  rounding  errors.   The
  234.     difference  is  however  the  smallest  possible,  i.e.,  less than
  235.     BSCALE / 2.
  236.  
  237.     Keywords DATA_MIN and DATA_MAX may be used  to  supply  the maximum
  238.     and minimum data values or to set brightness cutoffs.
  239.  
  240.    SEE ALSO fitsRead, fitsWrite.
  241. */
  242. {
  243.   if (bitpix == 8) {
  244.     // assume ``char'' is 8 bits unsigned
  245.     file_type= char;
  246.     file_unsigned= 1N;
  247.   } else if (bitpix == 16) {
  248.     // assume ``short'' is 16 bits signed
  249.     file_type= short;
  250.     file_unsigned= 0N;
  251.   } else if (bitpix == 32) {
  252.     // assume ``long'' is 32 bits signed
  253.     file_type= long;
  254.     file_unsigned= 0N;
  255.   } else {
  256.     error, "bad BITPIX (should be 8, 16 or 32)";
  257.   }
  258.   data_min= double((is_void(data_min) ? min(data) : data_min));
  259.   data_max= double((is_void(data_max) ? max(data) : data_max));
  260.   if (data_max < data_min)
  261.     error, "bad DATA_MAX and DATA_MIN";
  262.   if (file_unsigned) {
  263.     file_min= 0.;
  264.     file_max= 2.^bitpix - 1.;
  265.   } else {
  266.     file_min= -2.^(bitpix - 1);
  267.     file_max=  2.^(bitpix - 1) - 1.;
  268.   }
  269.   if (data_max == data_min) {
  270.     bzero= data_min;
  271.     return array(file_type(0), dimsof(data));
  272.   }
  273.   bscale= (data_max - data_min) / (file_max - file_min);
  274.   bzero= data_min - bscale * file_min;
  275.   return file_type(min(file_max, max(file_min, (data - bzero) / bscale)));
  276. }
  277.  
  278. /*----------------------------------------------------------------------*/
  279.  
  280. func fitsWrite(name, data, header, rescale=, pack=)
  281. /* DOCUMENT fitsWrite, filename, data, header
  282.     Write DATA in file FILENAME using  FITS  format.   If  present, the
  283.     information  of  the  optional   argument   HEADER   (a  FitsHeader
  284.     structure) will be used to write the FITS file header.
  285.  
  286.     Keyword PACK, if non-nil and non-zero,  indicates  that  axis whith
  287.     unit dimension should be ignored.  The default  is  to  ignore only
  288.     zero length axis.
  289.  
  290.     Keyword RESCALE, if non-nil and  zero,  indicates  that  input data
  291.     values should not  be  rescaled  into  integers  according  to FITS
  292.     standard.  The default, is to recode doubles  as  32  bit integers,
  293.     floats as 16 bit integers and integers into smallest integers if it
  294.     is possible  without  loss  of  dynamic  (e.g.,  an  array  of long
  295.     integers ranging from -31 to +130 will be recoded as chars).
  296.  
  297.    SEE ALSO: fitsRead, fitsRescale, fitsAddHistory, fitsAddComment.
  298. */
  299. {
  300.   local bitpix, bzero, bscale;
  301.  
  302.   /*
  303.    * Check input data.
  304.    */
  305.   if (is_void(data))
  306.     error, "empty DATA array!!!";
  307.   if (dimsof(data)(1) > fits_max_naxis) {
  308.     write, swrite(format="*** Currently, number of dimensions for "+
  309.           "FITS cannot exceed %d\n*** If you really need "+
  310.           "more dimensions, you'll have to modify the\n"+
  311.           "*** source file...", int(fits_max_naxis));
  312.     error, "too many dimensions";
  313.   }
  314.   bitpix= fitsBitpix(data);    // also make sure that data type is valid
  315.  
  316.   /*
  317.    * Minimal header.
  318.    */
  319.   if (is_void(header)) {
  320.     fitsHeader, header;    // create a FITS header with defaults
  321.   } else {
  322.     fitsFixHeader, header;    // just check header
  323.   }
  324.   if (strlen(header.date) == 0) {
  325.     header.date= getdate()+" "+gettime();
  326.   }
  327.   header.bitpix= bitpix;
  328.   if (header.bscale != 1. || header.bzero != 0.)
  329.     write, "WARNING: initial [BSCALE,BZERO] != [1,0]";
  330.  
  331.   /*
  332.    * Rescale brightness to fit standards.
  333.    */
  334.   if (is_void(rescale) || rescale) {
  335.     if (structof(data) == float || structof(data) == double) {
  336.       // Convert float to short and double to long
  337.       bitpix= fitsBitpix((structof(data) == float ? short(0) : long(0)));
  338.       data= fitsRescale(data, bitpix, bscale, bzero);
  339.       header.bitpix = bitpix;
  340.       header.bzero += header.bscale * bzero;
  341.       header.bscale *= bscale;
  342.     } else if (structof(data) != char) {
  343.       // Integers: check if data fit into smaller data type
  344.       data_min= double(min(data));
  345.       data_max= double(max(data));
  346.       data_dyn= long(ceil(log(data_max - data_min + 1.) / log(256.)));
  347.       if (data_dyn < sizeof(data(1))) {
  348.     if (data_dyn == sizeof(char)) {
  349.       file_type= char;
  350.       file_unsigned= 1N;    // assume ``char'' is unsigned
  351.     } else {
  352.       file_type= short;
  353.       file_unsigned= 0N;    // assume ``short'' is signed
  354.     }
  355.     if (file_unsigned)
  356.       bzero= data_min;
  357.     else
  358.       bzero= data_min + 2.^(8. * data_siz - 1);
  359.     data= file_type(data - bzero);
  360.     header.bitpix= fitsBitpix(file_type(0));
  361.     header.bzero += header.bscale * bzero;
  362.       }
  363.     }
  364.   }
  365.  
  366.   /*
  367.    * Compute dimensions list.
  368.    */
  369.   header.axis(*)= 0;
  370.   dims= dimsof(data);
  371.   if (dims(1) == 0) {
  372.     // data is a scalar!
  373.     header.naxis= 1;
  374.     header.axis(1)= 1;
  375.   } else {
  376.     // there is at least one dimension
  377.     dims = dims(2:);
  378.     if (pack) {
  379.       i= where(dims > 1);
  380.       if (numberof(i)) {
  381.     dims= dims(i);
  382.       } else {
  383.     dims= [1];
  384.       }
  385.     }
  386.     header.naxis = numberof(dims);
  387.     header.axis(:header.naxis) = dims;
  388.   }
  389.  
  390.   /*
  391.    * Open data file.
  392.    */
  393.   if (open(name+"L", "", 1)) {
  394.     write, format="%s", "Remove existing file "+name+"L (y or n)? ";
  395.     response= strtok(rdline(prompt=string(0)))(1);
  396.     if (response=="y" || response=="yes") {
  397.       remove, name+"L";
  398.     } else {
  399.       write, "No file written or deleted.";
  400.       return;
  401.     }
  402.   }
  403.   file= open(name, "wb");
  404.   sun_primitives, file;    // Most Significant Byte First
  405.   address= 0;        // Address in file
  406.  
  407.  
  408.   /*
  409.    * Write header.
  410.    */
  411.   _fitsPutLogical, "SIMPLE", 1;
  412.   _fitsPutInteger, "BITPIX", header.bitpix;
  413.   _fitsPutInteger, "NAXIS", header.naxis;
  414.   for (i=1; i<=header.naxis; i++) {
  415.     _fitsPutInteger, swrite(format="NAXIS%d", i), header.axis(i);
  416.   }
  417.   if (header.bscale != 1.0 || header.bzero != 0.) {
  418.     _fitsPutReal, "BSCALE", header.bscale;
  419.     _fitsPutReal, "BZERO", header.bzero;
  420.   }
  421.   if (strlen(header.bunit))
  422.     _fitsPutString, "BUNIT", header.bunit;
  423.   if (header.datamax || header.datamin) {
  424.     _fitsPutReal, "DATAMAX", header.datamax;
  425.     _fitsPutReal, "DATAMIN", header.datamin;
  426.   }
  427.   if (strlen(header.object))
  428.     _fitsPutString, "OBJECT", header.object;
  429.   if (strlen(header.date))
  430.     _fitsPutString, "DATE", header.date;
  431.   if (strlen(header.date_obs))
  432.     _fitsPutString, "DATE-OBS", header.date_obs;
  433.   if (strlen(header.origin))
  434.     _fitsPutString, "ORIGIN", header.origin;
  435.   if (strlen(header.instrume))
  436.     _fitsPutString, "INSTRUME", header.instrume;
  437.   if (strlen(header.telescop))
  438.     _fitsPutString, "TELESCOP", header.telescop;
  439.   if (strlen(header.observer))
  440.     _fitsPutString, "OBSERVER", header.observer;
  441.   if (header.epoch)
  442.     _fitsPutReal, "EPOCH", header.epoch;
  443.   for (i=1; i<=header.naxis; i++) {
  444.     if (strlen(header.ctype(i)))
  445.       _fitsPutString, swrite(format="CTYPE%d", i), header.ctype(i);
  446.     if (header.cdelt(i)) {
  447.       _fitsPutReal, swrite(format="CRVAL%d", i), header.crval(i);
  448.       _fitsPutReal, swrite(format="CRPIX%d", i), header.crpix(i);
  449.       _fitsPutReal, swrite(format="CDELT%d", i), header.cdelt(i);
  450.       _fitsPutReal, swrite(format="CROTA%d", i), header.crota(i);
  451.     }
  452.   }
  453.   if (strlen(header.comment)) {
  454.     comment = strtok(header.comment, "\n");
  455.     while (comment(1)) {
  456.       _fitsPutComment, "COMMENT", comment(1);
  457.       comment = strtok(comment(2), "\n");
  458.     }
  459.   }
  460.   if (strlen(header.history)) {
  461.     history = strtok(header.history, "\n");
  462.     while (history(1)) {
  463.       _fitsPutComment, "HISTORY", history(1);
  464.       history = strtok(history(2), "\n");
  465.     }
  466.   }
  467.  
  468.   /*
  469.    * Write END keyword and terminate current header block.
  470.    */
  471.   _fitsPutComment, "END";
  472.   empty= array(' ', 80);
  473.   empty(0)= '\n';
  474.   while (address % 2880) {
  475.     _write, file, address, empty;
  476.     address += 80;
  477.   }
  478.  
  479.   /*
  480.    * Write data array.
  481.    */
  482.   _write, file, address, data;
  483.   close, file;
  484.  
  485.   /*
  486.    * Yorick uses a Contents Log file, which should be removed...
  487.    */
  488.   if (open(name+"L", "", 1))
  489.     remove, name+"L";
  490. }
  491.  
  492. /*----------------------------------------------------------------------*/
  493.  
  494. func fitsRead(name, &header, which=, pack=, rescale=)
  495. /* DOCUMENT a= fitsRead(filename, header)
  496.  
  497.     Returns the data  of  the  FITS  file  FILENAME.   If  present, the
  498.     optional argument HEADER will be used to store the contents  of the
  499.     FITS header file (a FitsHeader structure).
  500.  
  501.     Keyword WHICH may be used to  indicate  which  sub-array  should be
  502.     returned.   For  instance,  if  the  array  DATA   with  dimensions
  503.     (235,453,7) is stored in the FITS file  "data.fits",  the sub-array
  504.     DATA(,,4) can be read by:
  505.         SUB_DATA= fitsRead("data.fits", which= 4);
  506.  
  507.     Keyword PACK, if non-nil and non-zero,  indicates  that  axis whith
  508.     unit dimension should be ignored.  The default  is  to  ignore only
  509.     zero length axis.
  510.  
  511.     Keyword RESCALE, if non-nil  and  zero,  indicates  that  read data
  512.     values should not be rescaled according to FITS keywords BSCALE and
  513.     BZERO.  The default is to rescale  data  values  if  BSCALE  is not
  514.     1. or BZERO is not 0.
  515.  
  516.    SEE ALSO: fitsWrite, fitsRescale.
  517. */
  518. {
  519.   local axis;
  520.  
  521.   fitsHeader, header;        // Create FITS header with defaults
  522.   file= open(name, "rb");
  523.   sun_primitives, file;    // Most Signficant Byte First
  524.   end_not_found= 1N;
  525.   n= 1;
  526.   address= 0;
  527.   do {
  528.     /*
  529.      * Read header block.
  530.      */
  531.     buffer= array(char, 80, 36);
  532.     if (_read(file, address, buffer) != sizeof(buffer))
  533.       error, "cannot read header";
  534.     address += sizeof(buffer);
  535.     
  536.     /*
  537.      * Get keywords: convert to upper case letters and remove trailing
  538.      * blanks.
  539.      */
  540.     k= buffer(1:8,);
  541.     v= buffer(9:,);
  542.     if (numberof((i=where(k >= 'a'))) && numberof((j=where(k(i) <= 'z'))))
  543.       k(i(j))+= long('A')-long('a');
  544.     keyword= array(string, 36);
  545.     value= array(string, 36);
  546.     for (i=1; i<=36; i++) {
  547.       sread, string(&k(,i)), format="%[^ \t]", keyword(i);
  548.       value(i)= string(&v(,i));
  549.     }
  550.     k= v= buffer= [];
  551.  
  552.     /*
  553.      * Parse minimal header.
  554.      */
  555.     if (n > 1) {
  556.       i= 0;
  557.     } else {
  558.       if (keyword(1) != "SIMPLE" || !_fitsGetLogical(keyword(1), value(1)))
  559.     error, "not a standard FITS file";
  560.       if (keyword(2) != "BITPIX" ||
  561.       noneof(((header.bitpix= _fitsGetInteger(keyword(2), value(2))) ==
  562.           [8, 16, 32, 64, -32, -64, -128])))
  563.     error, "bad BITPIX or no BITPIX at line 2";
  564.       if (keyword(3) != "NAXIS" ||
  565.       (header.naxis= _fitsGetInteger(keyword(3), value(3))) < 1)
  566.     error, "bad NAXIS or no NAXIS at line 3";
  567.       if (header.naxis > fits_max_naxis) {
  568.     write, swrite(format="*** Currently, number of dimensions for "+
  569.               "FITS cannot exceed %d\n*** If you really need "+
  570.               "more dimensions, you'll have to modify the\n"+
  571.               "*** source file...", int(fits_max_naxis));
  572.     error, swrite(format="NAXIS=%d too large", header.naxis);
  573.       }
  574.       naxis= header.naxis;
  575.       for (k=1; k<=naxis; k++) {
  576.     i= k+3;
  577.     s= swrite(format="NAXIS%d", k);
  578.     if (keyword(i) != s ||
  579.         (header.axis(k)= _fitsGetInteger(keyword(i), value(i))) < 0)
  580.       error, "bad "+s;
  581.       }
  582.       n= i;
  583.     }
  584.     
  585.     /*
  586.      * Parse other keywords.
  587.      */
  588.     while (++i <= 36) {
  589.       n++;
  590.       k= keyword(i);
  591.       v= value(i);
  592.       if (strlen(k) == 0)
  593.     continue;
  594.       if (k == "END") {
  595.     end_not_found= 0N;
  596.     break;
  597.       }
  598.       if (k == "COMMENT") {
  599.     v= strtrim(v);
  600.     if (strlen(v))
  601.       header.comment += v+"\n";
  602.     continue;
  603.       }
  604.       if (k == "HISTORY") {
  605.     v= strtrim(v);
  606.     if (strlen(v))
  607.       header.history += v+"\n";
  608.     continue;
  609.       }
  610.       if (k == "BSCALE") {
  611.     header.bscale = _fitsGetReal(k, v);
  612.     continue;
  613.       }
  614.       if (k == "BZERO") {
  615.     header.bzero = _fitsGetReal(k, v);
  616.     continue;
  617.       }
  618.       if (k == "BUNIT") {
  619.     header.bunit = _fitsGetString(k, v);
  620.     continue;
  621.       }
  622.       if (k == "DATAMAX") {
  623.     header.datamax = _fitsGetReal(k, v);
  624.     continue;
  625.       }
  626.       if (k == "DATAMIN") {
  627.     header.datamin = _fitsGetReal(k, v);
  628.     continue;
  629.       }
  630.       if (k == "EPOCH") {
  631.     header.epoch = _fitsGetReal(k, v);
  632.     continue;
  633.       }
  634.       if (k == "OBJECT") {
  635.     header.object = _fitsGetString(k, v);
  636.     continue;
  637.       }
  638.       if (k == "DATE") {
  639.     header.date = _fitsGetString(k, v);
  640.     continue;
  641.       }
  642.       if (k == "DATE-OBS") {
  643.     header.date_obs = _fitsGetString(k, v);
  644.     continue;
  645.       }
  646.       if (k == "ORIGIN") {
  647.     header.origin = _fitsGetString(k, v);
  648.     continue;
  649.       }
  650.       if (k == "INSTRUME") {
  651.     header.instrume = _fitsGetString(k, v);
  652.     continue;
  653.       }
  654.       if (k == "TELESCOP") {
  655.     header.telescop = _fitsGetString(k, v);
  656.     continue;
  657.       }
  658.       if (k == "OBSERVER") {
  659.     header.observer = _fitsGetString(k, v);
  660.     continue;
  661.       }
  662.       if (k == "BLANK") {
  663.     if (strpart(v, 1:1) == "=") {
  664.       header.blank = _fitsGetInteger(k, v);
  665.     } else {
  666.       v= strtrim(v);
  667.       if (strlen(v))
  668.         write, "WARNING: discarding BLANK comment \""+v+"\"";
  669.     }
  670.     continue;
  671.       }
  672.       if (strpart(k, 1:1) == "C") {
  673.     j= long();
  674.     s= string();
  675.     if (sread(k, format="%[^0-9]%d", s, j) != 2)
  676.       continue;
  677.     if (s == "CRVAL") {
  678.       if (j < 1 || j > header.naxis)
  679.         error, "bad "+k;
  680.       header.crval(j)= _fitsGetReal(k, v);
  681.       continue;
  682.     }
  683.     if (s == "CRPIX") {
  684.       if (j < 1 || j > header.naxis)
  685.         error, "bad "+k;
  686.       header.crpix(j)= _fitsGetReal(k, v);
  687.       continue;
  688.     }
  689.     if (s == "CDELT") {
  690.       if (j < 1 || j > header.naxis)
  691.         error, "bad "+k;
  692.       header.cdelt(j)= _fitsGetReal(k, v);
  693.       continue;
  694.     }
  695.     if (s == "CROTA") {
  696.       if (j < 1 || j > header.naxis)
  697.         error, "bad "+k;
  698.       header.crota(j)= _fitsGetReal(k, v);
  699.       continue;
  700.     }
  701.     if (s == "CTYPE") {
  702.       if (j < 1 || j > header.naxis)
  703.         error, "bad "+k;
  704.       header.ctype(j)= _fitsGetString(k, v);
  705.       continue;
  706.     }
  707.       }
  708.       write, "WARNING: unknown keyword \""+k+"\"";
  709.     }
  710.   } while (end_not_found);
  711.   
  712.   /*
  713.    * Get dimensions (ignore zero and, possibly, unit dimensions, keep
  714.    * unit dimensions if keyword `pack' is zero or nil):
  715.    */
  716.   i= where(header.axis(:header.naxis) > (pack ? 1 : 0));
  717.   dims=  numberof(i);
  718.   if (dims <= 0) {
  719.     write, "WARNING empty data file";
  720.     close, file;
  721.     return [];
  722.   }
  723.   grow, dims, header.axis(i);
  724.  
  725.   /*
  726.    * Parse WHICH keyword for sub-array.
  727.    */
  728.   if (is_void(which)) {
  729.     which=0;
  730.   } else {
  731.     if (numberof(which) != 1 || which != long(which))
  732.       error, "WHICH must be a scalar integer";
  733.     last= dims(0);
  734.     if (which <= 0)
  735.       which += last;
  736.     if (which > last || which < 1)
  737.       error, "WHICH out of range";
  738.     dims= dims(:-1);
  739.     dims(1) -= 1;
  740.   }
  741.  
  742.   if (header.bitpix == 8) {
  743.     data = array(char, dims);
  744.   } else if (header.bitpix == 16) {
  745.     data = array(short, dims);
  746.   } else if (header.bitpix == 32) {
  747.     data = array(long, dims);
  748.   } else if (header.bitpix == -32) {
  749.     data = array(float, dims);
  750.   } else if (header.bitpix == -64) {
  751.     data = array(double, dims);
  752.   } else if (header.bitpix == -128) {
  753.     write, "WARNING: using BITPIX=-128 is not standard FITS";
  754.     data = array(complex, dims);
  755.   } else {
  756.     error, "congratulation: you have evidenced a BUG!";
  757.   }
  758.  
  759.   if (which > 1)
  760.     address += (which - 1) * sizeof(data);
  761.   if (_read(file, address, data) != numberof(data))
  762.     error, "cannot read data";
  763.   close, file;
  764.   fitsFixHeader, header;
  765.  
  766.  
  767.   /*
  768.    * Maybe rescale pixel values
  769.    */
  770.   if ((is_void(rescale) || rescale) &&
  771.       (header.bscale != 1. || header.bzero != 0.)) {
  772.     if (abs(header.bitpix) < 32) {
  773.       data=  float(header.bzero) + float(header.bscale) * float(data);
  774.     } else {
  775.       data= header.bzero + header.bscale * data;
  776.     }
  777.     header.bitpix= fitsBitpix(data);
  778.     header.bscale= 1.;
  779.     header.bzero= 0.;
  780.   }
  781.  
  782.   return data;
  783. }
  784.  
  785. /*----------------------------------------------------------------------*/
  786.  
  787. func _fitsGetLogical(keyword, value)
  788. {
  789.   x = string();
  790.   if (strpart(value, 1:1) == "=" &&
  791.       sread(value, format="=%s", x) == 1 && strlen(x) == 1) {
  792.     if (x == "T" || x == "t")
  793.       return 1;
  794.     if (x == "F" || x == "f")
  795.       return 0;
  796.   }
  797.   error, "bad logical value for "+keyword;
  798. }
  799.  
  800. func _fitsGetInteger(keyword, value)
  801. {
  802.   x = long();
  803.   if (strpart(value, 1:1) == "=" && sread(value, format="=%d", x) == 1)
  804.     return x;
  805.   error, "bad integer value for "+keyword;
  806. }
  807.  
  808. func _fitsGetReal(keyword, value)
  809. {
  810.   x = double();
  811.   if (strpart(value, 1:1) == "=" && sread(value, format="=%f", x) == 1)
  812.     return x;
  813.   error, "bad real value for "+keyword;
  814. }
  815.  
  816. func _fitsGetString(keyword, value)
  817. {
  818.   if (strpart(value, 1:1) == "=") {
  819.     /* note that this does not allow for a / delimited comment field... */
  820.     x= strtrim(strpart(value, 2:));
  821.     if (strpart(x, 1:1) == "'") {
  822.       t= string();
  823.       sread, strpart(x, 2:), format="%[^']", t;
  824.       return strtrim(t);
  825.     } else {
  826.       return x;
  827.     }
  828.   }
  829.   //error, "bad value for "+keyword;
  830.   write, "WARNING: bad string value for "+keyword;
  831.   return "";
  832. }
  833.  
  834. /*----------------------------------------------------------------------*/
  835.  
  836. func _fitsPutValue(fmt, keyword, value, comment)
  837. {
  838.   extern file, address;
  839.   comment= is_void(comment) ? "" : strpart(comment, :46);
  840.   line= swrite(format=fmt, keyword, value, comment);
  841.   if (strlen(line) != 80)
  842.     error, "(BUG) bad line length: \""+line+"\"";
  843.   _write, file, address, (*pointer(line))(:80);
  844.   address += 80;
  845. }
  846.  
  847. func _fitsPutComment(keyword, comment)
  848. {
  849.   extern file, address;
  850.   comment= is_void(comment) ? "" : strpart(comment, 1:70);
  851.   line= swrite(format="%-8s %-70s\n", keyword, comment);
  852.   _write, file, address, (*pointer(line))(:80);
  853.   address += 80;
  854. }
  855.  
  856. func _fitsPutInteger(keyword, value, comment)
  857. {
  858.   _fitsPutValue, "%-8s= %20d / %-46s\n", keyword, long(value), comment;
  859. }
  860.  
  861. func _fitsPutReal(keyword, value, comment)
  862. {
  863.   _fitsPutValue, "%-8s= %#20.13G / %-46s\n", keyword, double(value), comment;
  864. }
  865.  
  866. func _fitsPutString(keyword, value, comment)
  867. {
  868.   _fitsPutValue, "%-8s= '%-18s' / %-46s\n", keyword, strpart(value, :18),
  869.     comment;
  870. }
  871.  
  872. func _fitsPutLogical(keyword, value, comment)
  873. {
  874.   _fitsPutValue, "%-8s= %20s / %-46s\n", keyword, (value ? "T" : "F"), comment;
  875. }
  876.  
  877. /*----------------------------------------------------------------------*/
  878.